#############################################################################################
########### calculate exposure inex for each species ########################################
#############################################################################################

#load packages:
library(R.utils)
library(raster)
library(sp)
library(sf)
library(stringr)
library(circular)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
MUOutDir<-paste(MyDir, "/Outputs2_Crossval/permutation.importance", sep="")
IndSppOutDir<-paste(MyDir, "/Outputs3_Final/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define 'for' variables:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
timecats<-c("2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
#timecats<-c("Current","2050_Median","2090_Median")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-DatSum[DatSum$Listing_comment!="not listed",]
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]
spSTART<-which (spp=="Vipera_ammodytes")

DatSumMod<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
# DatSumMod<-DatSumMod[,c("SSN","Family","Accepted_Species","Listing_Name","ModUnit","MergeSpecies","MergeShort","country")]
DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]



# # make extra columns to fill
# DatSumMod$SPrecs <- NA
# DatSumMod$MUrecs <- NA
# DatSumMod$TrainingAUC <- NA
# DatSumMod$TestAUC <- NA
# 
# DatSumMod$RangesizeSQKM_Current <- NA
# DatSumMod$SuitabilitySum_Current <- NA
# 
# DatSumMod$RangesizeSQKM_2050_Median <- NA
# DatSumMod$SuitabilitySum_2050_Median <- NA
# DatSumMod$RangesizeSQKM_2090_Median <- NA
# DatSumMod$SuitabilitySum_2090_Median <- NA
# 
# DatSumMod$RangesizeSQKM_2050_Q10 <- NA
# DatSumMod$SuitabilitySum_2050_Q10 <- NA
# DatSumMod$RangesizeSQKM_2090_Q10 <- NA
# DatSumMod$SuitabilitySum_2090_Q10 <- NA
# 
# DatSumMod$RangesizeSQKM_2050_Q90 <- NA
# DatSumMod$SuitabilitySum_2050_Q90 <- NA
# DatSumMod$RangesizeSQKM_2090_Q90 <- NA
# DatSumMod$SuitabilitySum_2090_Q90 <- NA
# 
# DatSumMod$lat_Current <- NA
# DatSumMod$lat_2050_Median <- NA
# DatSumMod$lat_2090_Median <- NA
# DatSumMod$lat_2050_Q10 <- NA
# DatSumMod$lat_2050_Q90 <- NA
# DatSumMod$lat_2090_Q10 <- NA
# DatSumMod$lat_2090_Q90 <- NA
# 
# DatSumMod$long_Current <- NA
# DatSumMod$long_2050_Median <- NA
# DatSumMod$long_2090_Median <- NA
# DatSumMod$long_2050_Q10 <- NA
# DatSumMod$long_2050_Q90 <- NA
# DatSumMod$long_2090_Q10 <- NA
# DatSumMod$long_2090_Q90 <- NA
# 
# DatSumMod$Risk_Current <- NA
# DatSumMod$Exposure_Current <- NA
# DatSumMod$Impact_Current <- NA
# 
# DatSumMod$Risk_2050_Median <- NA
# DatSumMod$Exposure_2050_Median <- NA
# DatSumMod$Impact_2050_Median <- NA
# DatSumMod$Risk_2090_Median <- NA
# DatSumMod$Exposure_2090_Median <- NA
# DatSumMod$Impact_2090_Median <- NA
# 
# DatSumMod$Risk_2050_Q10 <- NA
# DatSumMod$Exposure_2050_Q10 <- NA
# DatSumMod$Impact_2050_Q10 <- NA
# DatSumMod$Risk_2090_Q10 <- NA
# DatSumMod$Exposure_2090_Q10 <- NA
# DatSumMod$Impact_2090_Q10 <- NA
# 
# DatSumMod$Risk_2050_Q90 <- NA
# DatSumMod$Exposure_2050_Q90 <- NA
# DatSumMod$Impact_2050_Q90 <- NA
# DatSumMod$Risk_2090_Q90 <- NA
# DatSumMod$Exposure_2090_Q90 <- NA
# DatSumMod$Impact_2090_Q90 <- NA
# 
# 
# 
# DatSumMod$Loss_2050_Median <- NA
# DatSumMod$Loss_2090_Median <- NA
# DatSumMod$Loss_2050_Q10 <- NA
# DatSumMod$Loss_2050_Q90 <- NA
# DatSumMod$Loss_2090_Q10 <- NA
# DatSumMod$Loss_2090_Q90 <- NA
# 
# DatSumMod$Gain_2050_Median <- NA
# DatSumMod$Gain_2090_Median <- NA
# DatSumMod$Gain_2050_Q10 <- NA
# DatSumMod$Gain_2050_Q90 <- NA
# DatSumMod$Gain_2090_Q10 <- NA
# DatSumMod$Gain_2090_Q90 <- NA

# # read LN people per cell & % agriculture
# if(file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc.gz")) &
#    !file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"))){
#   gunzip(filename=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc.gz"),
#          destname=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"),
#          remove=FALSE, overwrite=TRUE)
# }
# if(file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc.gz")) &
#    !file.exists(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"))){
#   gunzip(filename=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc.gz"),
#          destname=paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"),
#          remove=FALSE, overwrite=TRUE)
# }
# 
# ppc<-raster(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/wpppcln.asc"), crs="+init=epsg:4326")
# pag<-raster(paste("/home/jc217070/GIS_and_raw_Layers/final_layers/percag.asc"), crs="+init=epsg:4326")

worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))


#for (sp in spp){
for (sp in spp[spSTART:length(spp)]){
  
  print(paste(Sys.time(), "appending table for", sp, sep=" "))

  print(paste(Sys.time(), "model stats", sep=" "))
  
  n<-which(spp==sp) #gets row number for current species
  mapsp<-DatSum$ModUnit[n]
  s<-which(DatSumMod$Accepted_Species==gsub("_"," ",sp)) #gets row number for current species
  
  
  
  
  
  
  if(file.exists(paste(OccDir, "/", sp,".csv",sep=""))){
  sprecs<-read.table(file=file.path(paste(OccDir, "/", sp,".csv",sep="")), header=TRUE, sep=",")
  NSPrecs<-length(sprecs[,1])
  # DatSumMod$SPrecs [s] <- NSPrecs
  }
  
  maprecs<-read.table(file=file.path(paste(OccDir, "/ModUnit_", gsub(" ","_",mapsp),".csv",sep="")), header=TRUE, sep=",")
  NMUrecs<-length(maprecs[,1])
  # DatSumMod$MUrecs [s] <- NMUrecs
  
  results <-read.csv(file=paste(MUOutDir, "/ModUnit_", gsub(" ","_",mapsp),"/maxentResults.csv", sep=""), header=TRUE)
  resultsB <-read.csv(file=paste(MyDir, "/Outputs3_Final/permutation.importance/ModUnit_", 
                                 gsub(" ","_",mapsp),"/maxentResults.csv", sep=""), header=TRUE)

  DatSumMod$TrainingAUC [s] <- results$Training.AUC[length(results$Training.AUC)]
  DatSumMod$TestAUC [s] <- results$Test.AUC[length(results$Test.AUC)]
  DatSumMod$TrainingAUC_FIN [s] <- resultsB$Training.AUC[length(resultsB$Training.AUC)]
  
# calculate other ?!!!stats!!!?, 

for (cat in timecats){
  
  print(paste(Sys.time(), "SDM stats", sep=" "))
  
#  outfile<-paste(IndSppOutDir,"/",sp,"_",cat, "_CropThreshCDCut.asc",sep="")
  
#   # unzip raster
#   if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
#   gunzip(filename=paste(outfile,".gz", sep=""),
#          destname=paste(outfile),
#          remove=FALSE, overwrite=TRUE)
#   }
#   
#   # read raster
#   SDM<-raster(paste(outfile), crs="+init=epsg:4326")
#   
#   
#   # get all values
#   myvals <- values(SDM)
#   myvals<- myvals[myvals>0]
#   myvals<- myvals[!is.na(myvals)]
#   mysum<-sum(myvals)
#   myrange <- length(myvals)
#   
#   
#   mycoords<-as.data.frame(xyFromCell(SDM, which(SDM[]>0)))
#   myX<-mean(mycoords$x)
#   myY<-mean(mycoords$y)
#   
#   # append table with sum (total suitability) & range size
#   mycolname<-paste("RangesizeSQKM_",cat,sep="")
#   myc<-which(names(DatSumMod)==mycolname)
#   DatSumMod [n,myc] <- myrange
#   
#   mycolname<-paste("SuitabilitySum_",cat,sep="")
#   myc<-which(names(DatSumMod)==mycolname)
#   DatSumMod [n,myc] <- mysum
#   
#   mycolname<-paste("lat_",cat,sep="")
#   myc<-which(names(DatSumMod)==mycolname)
#   DatSumMod [n,myc] <- myY
#   
#   mycolname<-paste("long_",cat,sep="")
#   myc<-which(names(DatSumMod)==mycolname)
#   DatSumMod [n,myc] <- myX
#   
#   ##############################################
#   
#   print(paste(Sys.time(), "Risk Index", sep=" "))
#   
# #Risk Index
# Risk<-pag*SDM
# 
# #SAVE!!!!#
# writeRaster(Risk, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Risk.asc",sep=""), format="ascii", overwrite=TRUE)
# gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Risk.asc",sep=""),overwrite=TRUE, remove=TRUE)
# 
# # get all values & sum
# myvals <- values(Risk)
# myvals<- myvals[myvals>0]
# myvals<- myvals[!is.na(myvals)]
# mysum<-sum(myvals)
# # append table row with sum (risk index)
# mycolname<-paste("Risk_",cat,sep="")
# myc<-which(names(DatSumMod)==mycolname)
# DatSumMod [n,myc] <- mysum
# 
# print(paste(Sys.time(), "Exposure Index", sep=" "))
# 
# #Exposure Index
# Exp<-ppc*SDM
# 
# 
# #SAVE!!!!#
# writeRaster(Exp, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Exp.asc",sep=""), format="ascii", overwrite=TRUE)
# gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Exp.asc",sep=""),overwrite=TRUE, remove=TRUE)
# 
# 
# myvals <- values(Exp)
# myvals<- myvals[myvals>0]
# myvals<- myvals[!is.na(myvals)]
# mysum<-sum(myvals)
# DatSumMod [n,myc+1] <- mysum
# 
# print(paste(Sys.time(), "Impact Index", sep=" "))
# 
# #Impact Index
# Imp<-Exp*pag
# 
# 
# #SAVE!!!!#
# writeRaster(Imp, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Imp.asc",sep=""), format="ascii", overwrite=TRUE)
# gzip(filename=paste(IndSppOutDir,"/",sp,"_",cat,"_Imp.asc",sep=""),overwrite=TRUE, remove=TRUE)
# 
# 
# myvals <- values(Imp)
# myvals<- myvals[myvals>0]
# myvals<- myvals[!is.na(myvals)]
# mysum<-sum(myvals)
# DatSumMod [n,myc+2] <- mysum
# 
# 
# 
# #remove unzipped raster or zip
# if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#   file.remove(outfile)
# }
# if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#   gzip(filename=paste(outfile))
# }

  {  #####################################################################
    # #read lost vs gained files for calculations
    # outfile<-paste(IndSppOutDir,"/",sp,"_","2050_Median_LossVsGain.asc",sep="")
    # # unzip raster
    # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    #   gunzip(filename=paste(outfile,".gz", sep=""),
    #          destname=paste(outfile),
    #          remove=FALSE, overwrite=TRUE)
    # }
    # LvG2050<-raster(paste(outfile), crs="+init=epsg:4326")
    # 
    # outfile<-paste(IndSppOutDir,"/",sp,"_","2090_Median_LossVsGain.asc",sep="")
    # # unzip raster
    # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    #   gunzip(filename=paste(outfile,".gz", sep=""),
    #          destname=paste(outfile),
    #          remove=FALSE, overwrite=TRUE)
    # }
    # LvG2090<-raster(paste(outfile), crs="+init=epsg:4326")
  }
  
  
  {
    
    
    
    # # calculate suitability lost & gained
    # outfile<-paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep="")
    # # unzip raster
    # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    #   gunzip(filename=paste(outfile,".gz", sep=""),
    #          destname=paste(outfile),
    #          remove=FALSE, overwrite=TRUE)
    # }
    # Suit1985<-raster(paste(outfile), crs="+init=epsg:4326")
    # 
    # outfile<-paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep="")
    # # unzip raster
    # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    #   gunzip(filename=paste(outfile,".gz", sep=""),
    #          destname=paste(outfile),
    #          remove=FALSE, overwrite=TRUE)
    # }
    # Suit2050<-raster(paste(outfile), crs="+init=epsg:4326")
    # 
    # outfile<-paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep="")
    # # unzip raster
    # if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    #   gunzip(filename=paste(outfile,".gz", sep=""),
    #          destname=paste(outfile),
    #          remove=FALSE, overwrite=TRUE)
    # }
    # Suit2090<-raster(paste(outfile), crs="+init=epsg:4326")
    # 
    # myras<-LvG2050*Suit2050 #(anything>0)
    # myvals <- values(myras)
    # myvals<- myvals[myvals>0]
    # myvals<- myvals[!is.na(myvals)]
    # suitgain2050<-sum(myvals)
    # 
    # myras<-LvG2090*Suit2090 #(anything>0)
    # myvals <- values(myras)
    # myvals<- myvals[myvals>0]
    # myvals<- myvals[!is.na(myvals)]
    # suitgain2090<-sum(myvals)
    # 
    # myras<-LvG2050*Suit1985 #(anything<0)
    # myvals <- values(myras)
    # myvals<- myvals[myvals>0]
    # myvals<- myvals[!is.na(myvals)]
    # suitloss2050<-sum(myvals)
    # 
    # myras<-LvG2090*Suit1985 #(anything<0)
    # myvals <- values(myras)
    # myvals<- myvals[myvals>0]
    # myvals<- myvals[!is.na(myvals)]
    # suitloss2090<-sum(myvals)
    # 
    # SLossMed<-c(0,
    #             suitloss2050/DatSumMod$SuitabilitySum_Current[n],
    #             suitloss2090/DatSumMod$SuitabilitySum_Current[n])
    # SLossMed[SLossMed==0]<-0.01 #avoid dividing by 0 for % change graphs
    # 
    # if(file.exists(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc.gz",sep="")) & 
    #    !file.exists(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep=""))){
    #   file.remove(paste(IndSppOutDir,"/",sp,"_","2050_Median_CropThreshCDCut.asc",sep=""))
    # }
    # if(file.exists(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc.gz",sep="")) & 
    #    !file.exists(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep=""))){
    #   file.remove(paste(IndSppOutDir,"/",sp,"_","2090_Median_CropThreshCDCut.asc",sep=""))
    # }
    # if(file.exists(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc.gz",sep="")) & 
    #    !file.exists(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep=""))){
    #   file.remove(paste(IndSppOutDir,"/",sp,"_","Current_CropThreshCDCut.asc",sep=""))
    # }
    # 
    
    
  }
  
  
  # calculate total increases and decreases in suitability? (sum of all positive and negative values in suitchange maps)
  
  
  #################################
  #################################
  #read in range loss/gain raster:
  
    outfile<-paste(IndSppOutDir,"/",sp,"_",cat, "_LossVsGain.asc",sep="")
  
    # unzip raster
    if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
    gunzip(filename=paste(outfile,".gz", sep=""),
           destname=paste(outfile),
           remove=FALSE, overwrite=TRUE)
    }

    # read raster
    SDM<-raster(paste(outfile), crs="+init=epsg:4326")

    # get all values
    myvals <- values(SDM)
    
    #calculate gained range 
    myvals<- myvals[myvals>0]
    myvals<- myvals[!is.na(myvals)]
    mygain <- length(myvals)

    # get all values
    myvals <- values(SDM)
    
    #calculate lost range
    myvals<- myvals[myvals<0]
    myvals<- myvals[!is.na(myvals)]
    myloss <- length(myvals)
    

    #add to summary file
     mycolname<-paste("Gain_",cat,sep="")
     myc<-which(names(DatSumMod)==mycolname)
     DatSumMod [n,myc] <- mygain
    
     mycolname<-paste("Loss_",cat,sep="")
     myc<-which(names(DatSumMod)==mycolname)
     DatSumMod [n,myc] <- myloss
     
    #remove unzipped raster or zip
    if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
      file.remove(outfile)
    }
    if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
      gzip(filename=paste(outfile))
    }

     #################################
     #################################
     
  
#write copy of updated table
write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""),
            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
}}



DatSumMod$RangesizeSQKM_CHANGE_2050_Median  <-  DatSumMod$RangesizeSQKM_2050_Median - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2050_Q10     <-  DatSumMod$RangesizeSQKM_2050_Q10 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2050_Q90     <-  DatSumMod$RangesizeSQKM_2050_Q90 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Median  <-  DatSumMod$RangesizeSQKM_2090_Median - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Q10     <-  DatSumMod$RangesizeSQKM_2090_Q10 - DatSumMod$RangesizeSQKM_Current
DatSumMod$RangesizeSQKM_CHANGE_2090_Q90     <- DatSumMod$RangesizeSQKM_2090_Q90 - DatSumMod$RangesizeSQKM_Current

DatSumMod$SuitabilitySum_CHANGE_2050_Median <- DatSumMod$SuitabilitySum_2050_Median - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2050_Q10    <- DatSumMod$SuitabilitySum_2050_Q10 - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2050_Q90    <- DatSumMod$SuitabilitySum_2050_Q90 - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Median <- DatSumMod$SuitabilitySum_2090_Median - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Q10    <- DatSumMod$SuitabilitySum_2090_Q10 - DatSumMod$SuitabilitySum_Current
DatSumMod$SuitabilitySum_CHANGE_2090_Q90    <- DatSumMod$SuitabilitySum_2090_Q90 - DatSumMod$SuitabilitySum_Current

DatSumMod$Risk_CHANGE_2050_Median           <-  DatSumMod$Risk_2050_Median -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2050_Q10              <-  DatSumMod$Risk_2050_Q10 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2050_Q90              <-  DatSumMod$Risk_2050_Q90 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Median           <- DatSumMod$Risk_2090_Median -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Q10              <- DatSumMod$Risk_2090_Q10 -  DatSumMod$Risk_Current
DatSumMod$Risk_CHANGE_2090_Q90              <- DatSumMod$Risk_2090_Q90 -  DatSumMod$Risk_Current

DatSumMod$Exposure_CHANGE_2050_Median       <- DatSumMod$Exposure_2050_Median -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2050_Q10          <- DatSumMod$Exposure_2050_Q10 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2050_Q90          <- DatSumMod$Exposure_2050_Q90 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Median       <- DatSumMod$Exposure_2090_Median -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Q10          <- DatSumMod$Exposure_2090_Q10 -  DatSumMod$Exposure_Current
DatSumMod$Exposure_CHANGE_2090_Q90          <- DatSumMod$Exposure_2090_Q90 -  DatSumMod$Exposure_Current

DatSumMod$Impact_CHANGE_2050_Median         <- DatSumMod$Impact_2050_Median -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2050_Q10            <- DatSumMod$Impact_2050_Q10 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2050_Q90            <- DatSumMod$Impact_2050_Q90 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Median         <- DatSumMod$Impact_2090_Median -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Q10            <- DatSumMod$Impact_2090_Q10 -  DatSumMod$Impact_Current
DatSumMod$Impact_CHANGE_2090_Q90            <- DatSumMod$Impact_2090_Q90 -  DatSumMod$Impact_Current

#calculate change vectors (distance and direction)

cols1<-grep("lat",names(DatSumMod))
cols2<-grep("long",names(DatSumMod))
mycols<-sort(c(cols1,cols2))

for (col in   mycols){
  
  DatSumMod[,col]  <- ifelse(DatSumMod[,col]== 0,NA,DatSumMod[,col])
  
}



#set up variables:
R<- 6371e3
latpiCurrent = DatSumMod$lat_Current * pi/180
latpi2050Med = DatSumMod$lat_2050_Median * pi/180
latpi2090Med = DatSumMod$lat_2090_Median * pi/180
latpi2050Q10 = DatSumMod$lat_2050_Q10 * pi/180
latpi2050Q90 = DatSumMod$lat_2050_Q90 * pi/180
latpi2090Q10 = DatSumMod$lat_2090_Q10 * pi/180
latpi2090Q90 = DatSumMod$lat_2090_Q90 * pi/180
longpiCurrent = DatSumMod$long_Current * pi/180
longpi2050Med = DatSumMod$long_2050_Median * pi/180
longpi2090Med = DatSumMod$long_2090_Median * pi/180
longpi2050Q10 = DatSumMod$long_2050_Q10 * pi/180
longpi2050Q90 = DatSumMod$long_2050_Q90 * pi/180
longpi2090Q10 = DatSumMod$long_2090_Q10 * pi/180
longpi2090Q90 = DatSumMod$long_2090_Q90 * pi/180

#calculate distances:
DatSumMod$DistM_2050_Median <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Med-latpiCurrent)/2) * sin((latpi2050Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Med) * sin((longpi2050Med-longpiCurrent)/2) * sin((longpi2050Med-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Med-latpiCurrent)/2) * sin((latpi2050Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Med) *  sin((longpi2050Med-longpiCurrent)/2) * sin((longpi2050Med-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2050_Q10 <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Q10-latpiCurrent)/2) * sin((latpi2050Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q10) * sin((longpi2050Q10-longpiCurrent)/2) * sin((longpi2050Q10-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Q10-latpiCurrent)/2) * sin((latpi2050Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q10) *  sin((longpi2050Q10-longpiCurrent)/2) * sin((longpi2050Q10-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2050_Q90 <-  R * (2*atan2(sqrt     
                                             (sin((latpi2050Q90-latpiCurrent)/2) * sin((latpi2050Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q90) * sin((longpi2050Q90-longpiCurrent)/2) * sin((longpi2050Q90-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2050Q90-latpiCurrent)/2) * sin((latpi2050Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2050Q90) *  sin((longpi2050Q90-longpiCurrent)/2) * sin((longpi2050Q90-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2090_Median <-  R * (2*atan2(sqrt     
                                             (sin((latpi2090Med-latpiCurrent)/2) * sin((latpi2090Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Med) * sin((longpi2090Med-longpiCurrent)/2) * sin((longpi2090Med-longpiCurrent)/2))   
                                             , sqrt(1-  
                                                      (sin((latpi2090Med-latpiCurrent)/2) * sin((latpi2090Med-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Med) *  sin((longpi2090Med-longpiCurrent)/2) * sin((longpi2090Med-longpiCurrent)/2)) 
                                             )))

DatSumMod$DistM_2090_Q10 <-  R * (2*atan2(sqrt     
                                          (sin((latpi2090Q10-latpiCurrent)/2) * sin((latpi2090Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q10) * sin((longpi2090Q10-longpiCurrent)/2) * sin((longpi2090Q10-longpiCurrent)/2))   
                                          , sqrt(1-  
                                                   (sin((latpi2090Q10-latpiCurrent)/2) * sin((latpi2090Q10-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q10) *  sin((longpi2090Q10-longpiCurrent)/2) * sin((longpi2090Q10-longpiCurrent)/2)) 
                                          )))

DatSumMod$DistM_2090_Q90 <-  R * (2*atan2(sqrt     
                                          (sin((latpi2090Q90-latpiCurrent)/2) * sin((latpi2090Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q90) * sin((longpi2090Q90-longpiCurrent)/2) * sin((longpi2090Q90-longpiCurrent)/2))   
                                          , sqrt(1-  
                                                   (sin((latpi2090Q90-latpiCurrent)/2) * sin((latpi2090Q90-latpiCurrent)/2) +  cos(latpiCurrent) * cos(latpi2090Q90) *  sin((longpi2090Q90-longpiCurrent)/2) * sin((longpi2090Q90-longpiCurrent)/2)) 
                                          )))

#calculate directions:
DatSumMod$DirDeg_2050_Median = (180/pi) * atan2(sin(longpi2050Med-longpiCurrent) * cos(latpi2050Med) , cos(latpiCurrent)*sin(latpi2050Med) - sin(latpiCurrent)*cos(latpi2050Med)*cos(longpi2050Med-longpiCurrent))
DatSumMod$DirDeg_2050_Q10 = (180/pi) * atan2(sin(longpi2050Q10-longpiCurrent) * cos(latpi2050Q10) , cos(latpiCurrent)*sin(latpi2050Q10) - sin(latpiCurrent)*cos(latpi2050Q10)*cos(longpi2050Q10-longpiCurrent))
DatSumMod$DirDeg_2050_Q90 = (180/pi) * atan2(sin(longpi2050Q90-longpiCurrent) * cos(latpi2050Q90) , cos(latpiCurrent)*sin(latpi2050Q90) - sin(latpiCurrent)*cos(latpi2050Q90)*cos(longpi2050Q90-longpiCurrent))
DatSumMod$DirDeg_2090_Median = (180/pi) * atan2(sin(longpi2090Med-longpiCurrent) * cos(latpi2090Med) , cos(latpiCurrent)*sin(latpi2090Med) - sin(latpiCurrent)*cos(latpi2090Med)*cos(longpi2090Med-longpiCurrent))
DatSumMod$DirDeg_2090_Q10 = (180/pi) * atan2(sin(longpi2090Q10-longpiCurrent) * cos(latpi2090Q10) , cos(latpiCurrent)*sin(latpi2090Q10) - sin(latpiCurrent)*cos(latpi2090Q10)*cos(longpi2090Q10-longpiCurrent))
DatSumMod$DirDeg_2090_Q90 = (180/pi) * atan2(sin(longpi2090Q90-longpiCurrent) * cos(latpi2090Q90) , cos(latpiCurrent)*sin(latpi2090Q90) - sin(latpiCurrent)*cos(latpi2090Q90)*cos(longpi2090Q90-longpiCurrent))
# + = east, - = west
# 0 = north, 90 = east, -90 = west, 180/-180 =south





#!!!!!!!!!!!!!!!change from mean to median!!! (or do both)




# calculate summary stats for each value across all species for each time category

DatSumMod<-DatSumMod[DatSumMod$Accepted_Species!="Average",]
DatSumMod[length(DatSumMod[,1])+1,"Accepted_Species"]<- "Average"

for (col in   c(9:(length(DatSumMod)-6))  ){

  DatSumMod[length(DatSumMod[,1]),col]  <-mean(DatSumMod[1:(length(DatSumMod$Accepted_Species)-1),col],na.rm=TRUE)

  }

#do the circular mean for directions...
for (col in   c((length(DatSumMod)-5):length(DatSumMod))){
  #DatSumMod[length(DatSumMod[,1]),col]  <-mean(DatSumMod[,col],na.rm=TRUE)
  
  mydata  <-DatSumMod[1:(length(DatSumMod$Accepted_Species)-1),col]
  mydataC<-circular(mydata, type = "directions", 
                    units = "degrees",
                    modulo = "asis", 
                    rotation = "clock")
  
  DatSumMod[length(DatSumMod[,1]),col]<-mean.circular(mydataC,na.rm=TRUE)
}

write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""),
            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")



